library(lattice)
library(MASS)
library(e1071)
library(latticeExtra)
## Loading required package: RColorBrewer
library(corrplot)
concr.data <- read.csv("data/Concrete_Data.csv", comment.char = "#")
panel <- function(...) {
panel.xyplot(...)
panel.loess(...)
}
xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$Superplasticizer, panel=panel)
xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$Cement, panel=panel)
Видна отличная линейная зависимость прочности от цемента, а также от Superplasticizer(далее Суперпластик). Ну это и понятно, т.к. они являются главными скрепляющими и упрочняющими ингредиентами
xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$Water, panel=panel)
С водой происходит не очень понятный скачек вверх после систематического опускания. Вообще говоря мы достигаем минимума при самой большой плотности точек, а далее точек совсем мало, и возможно, дальшейшие результаты можно сбросить на недостаточное количество данных. Попробуем его позже разобрать подробнее.
xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$CoarseAggregate, panel=panel)
xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$FineAggregate, panel=panel)
Аггрегаты очень похожы, и переваливая за 810 практически не влияют на прочность, скорее всего они будут малозначимыми
xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$Age, panel=panel)
xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$BlastFurnaceSlag, panel=panel)
xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$FlyAsh, panel=panel)
Возраст выглядит совсем неоднородным, видимо его нужно будет сделать фактором
model1 <- lm(ConcreteCompressiveStrength ~ (.), data=concr.data)
summary(model1)
##
## Call:
## lm(formula = ConcreteCompressiveStrength ~ (.), data = concr.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.654 -6.302 0.703 6.569 34.450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -23.331214 26.585504 -0.878 0.380372
## Cement 0.119804 0.008489 14.113 < 2e-16 ***
## BlastFurnaceSlag 0.103866 0.010136 10.247 < 2e-16 ***
## FlyAsh 0.087934 0.012583 6.988 5.02e-12 ***
## Water -0.149918 0.040177 -3.731 0.000201 ***
## Superplasticizer 0.292225 0.093424 3.128 0.001810 **
## CoarseAggregate 0.018086 0.009392 1.926 0.054425 .
## FineAggregate 0.020190 0.010702 1.887 0.059491 .
## Age 0.114222 0.005427 21.046 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.4 on 1021 degrees of freedom
## Multiple R-squared: 0.6155, Adjusted R-squared: 0.6125
## F-statistic: 204.3 on 8 and 1021 DF, p-value: < 2.2e-16
Посмотрим на корреляцию:
corrplot(cor(concr.data))
xyplot(concr.data$Water ~ concr.data$Superplasticizer, panel=panel)
Построим модель, откинув не сильно значимые признаки агрегаторов, а также добавим Water:Superplasticizer, т.к. между ними явно есть зависимость
model1 <- lm(ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag + FlyAsh + Superplasticizer + Age + Water:Superplasticizer, data=concr.data)
summary(model1)
##
## Call:
## lm(formula = ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag +
## FlyAsh + Superplasticizer + Age + Water:Superplasticizer,
## data = concr.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.439 -6.400 0.209 6.940 33.095
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.502486 5.121626 8.884 < 2e-16 ***
## Cement 0.104717 0.004189 25.001 < 2e-16 ***
## Water -0.304965 0.026125 -11.673 < 2e-16 ***
## BlastFurnaceSlag 0.082116 0.004967 16.532 < 2e-16 ***
## FlyAsh 0.048647 0.008447 5.759 1.12e-08 ***
## Superplasticizer -2.351854 0.477719 -4.923 9.93e-07 ***
## Age 0.120990 0.005502 21.989 < 2e-16 ***
## Water:Superplasticizer 0.015927 0.002890 5.511 4.52e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.26 on 1022 degrees of freedom
## Multiple R-squared: 0.6252, Adjusted R-squared: 0.6226
## F-statistic: 243.5 on 7 and 1022 DF, p-value: < 2.2e-16
Немного лучше, но еще не очень.
Посмотрим на графики распределений предикторов:
marginal.plot(concr.data)
Age выглядит совсем плохо, попробуем прологарифмировать
marginal.plot(log(concr.data$Age))
xyplot(concr.data$ConcreteCompressiveStrength ~ log(concr.data$Age), panel=panel)
Стало лучше Посмотрим модель
model2 <- lm(ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag + FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer, data=concr.data)
summary(model2)
##
## Call:
## lm(formula = ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag +
## FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer,
## data = concr.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.9591 -4.4094 0.0202 4.2865 28.9493
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.464258 3.468202 5.612 2.57e-08 ***
## Cement 0.109381 0.002929 37.343 < 2e-16 ***
## Water -0.283281 0.017676 -16.027 < 2e-16 ***
## BlastFurnaceSlag 0.083428 0.003473 24.023 < 2e-16 ***
## FlyAsh 0.050869 0.005905 8.614 < 2e-16 ***
## Superplasticizer -1.164789 0.325079 -3.583 0.000356 ***
## log(Age) 8.712859 0.192143 45.346 < 2e-16 ***
## Water:Superplasticizer 0.007442 0.001965 3.787 0.000161 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.177 on 1022 degrees of freedom
## Multiple R-squared: 0.8167, Adjusted R-squared: 0.8154
## F-statistic: 650.4 on 7 and 1022 DF, p-value: < 2.2e-16
Намного лучше, но Мне все еще не нравится изломанность графика воды Поэтому можно попробовать добавить фактор для воды:
concr.data$Water.factor <- concr.data$Water < 225
model3 <-lm(ConcreteCompressiveStrength ~ Cement + Water.factor:Water + BlastFurnaceSlag + FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer, data=concr.data)
summary(model3)
##
## Call:
## lm(formula = ConcreteCompressiveStrength ~ Cement + Water.factor:Water +
## BlastFurnaceSlag + FlyAsh + Superplasticizer + log(Age) +
## Water:Superplasticizer, data = concr.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.7822 -4.3935 0.0625 4.5636 26.9560
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.895894 4.578324 8.059 2.14e-15 ***
## Cement 0.106611 0.002925 36.446 < 2e-16 ***
## BlastFurnaceSlag 0.080584 0.003456 23.316 < 2e-16 ***
## FlyAsh 0.042757 0.005986 7.143 1.74e-12 ***
## Superplasticizer -2.207334 0.368434 -5.991 2.89e-09 ***
## log(Age) 8.471278 0.193891 43.691 < 2e-16 ***
## Water.factorFALSE:Water -0.334456 0.019574 -17.087 < 2e-16 ***
## Water.factorTRUE:Water -0.368904 0.022963 -16.065 < 2e-16 ***
## Water:Superplasticizer 0.013515 0.002207 6.122 1.31e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.069 on 1021 degrees of freedom
## Multiple R-squared: 0.8224, Adjusted R-squared: 0.821
## F-statistic: 590.8 on 8 and 1021 DF, p-value: < 2.2e-16
Немного поперебирав значение порога получаем высоко значимый признак Water.factorTrue:Water, который говорит как влияет количство воды до 225 на прочность
Попробуем также и возраст сделать фактором:
concr.data$Age.factor <- concr.data$Age < 150
model4 <-lm(ConcreteCompressiveStrength ~ Cement + Water.factor:Water + BlastFurnaceSlag + FlyAsh + Superplasticizer + Age.factor:log(Age) + Water:Superplasticizer, data=concr.data)
summary(model4)
##
## Call:
## lm(formula = ConcreteCompressiveStrength ~ Cement + Water.factor:Water +
## BlastFurnaceSlag + FlyAsh + Superplasticizer + Age.factor:log(Age) +
## Water:Superplasticizer, data = concr.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.7246 -4.3742 -0.1258 4.3747 26.5782
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.336381 4.433372 7.294 6.04e-13 ***
## Cement 0.107538 0.002816 38.185 < 2e-16 ***
## BlastFurnaceSlag 0.078804 0.003331 23.657 < 2e-16 ***
## FlyAsh 0.040766 0.005764 7.073 2.82e-12 ***
## Superplasticizer -1.907606 0.356011 -5.358 1.04e-07 ***
## Water.factorFALSE:Water -0.306365 0.019084 -16.054 < 2e-16 ***
## Water.factorTRUE:Water -0.353940 0.022154 -15.976 < 2e-16 ***
## Age.factorFALSE:log(Age) 7.331764 0.224633 32.639 < 2e-16 ***
## Age.factorTRUE:log(Age) 9.301197 0.207623 44.799 < 2e-16 ***
## Water:Superplasticizer 0.011394 0.002137 5.333 1.19e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.801 on 1020 degrees of freedom
## Multiple R-squared: 0.8357, Adjusted R-squared: 0.8343
## F-statistic: 576.5 on 9 and 1020 DF, p-value: < 2.2e-16
И это тоже улучшило модель
Проверим модели
anova(model1, model2, model3, model4)
## Analysis of Variance Table
##
## Model 1: ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag +
## FlyAsh + Superplasticizer + Age + Water:Superplasticizer
## Model 2: ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag +
## FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer
## Model 3: ConcreteCompressiveStrength ~ Cement + Water.factor:Water + BlastFurnaceSlag +
## FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer
## Model 4: ConcreteCompressiveStrength ~ Cement + Water.factor:Water + BlastFurnaceSlag +
## FlyAsh + Superplasticizer + Age.factor:log(Age) + Water:Superplasticizer
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1022 107645
## 2 1022 52647 0 54998
## 3 1021 51013 1 1633 35.316 3.843e-09 ***
## 4 1020 47178 1 3835 82.924 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Последняя модель самая лучшая
xyplot(residuals(model1) ~ fitted(model1), panel=panel)
xyplot(residuals(model2) ~ fitted(model2), panel=panel)
xyplot(residuals(model3) ~ fitted(model3), panel=panel)
xyplot(residuals(model4) ~ fitted(model4), panel=panel)
Более менее линейно
cv <- function(model) {
tune(lm, model$call$formula, data = concr.data, tunecontrol = tune.control(sampling = "cross"))
}
cv(model1)
##
## Error estimation of 'lm' using 10-fold cross validation: 107.0302
cv(model2)
##
## Error estimation of 'lm' using 10-fold cross validation: 52.16315
cv(model3)
##
## Error estimation of 'lm' using 10-fold cross validation: 50.36774
cv(model4)
##
## Error estimation of 'lm' using 10-fold cross validation: 46.80454
На убывание, однако итоговая ошибка все-равно очень велика относительно диапозона прочности